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ABSTRACT 

We present 3D smoothed particle hydrodynamics simulations of the collapse of clumps 
formed through gravitational instability in the outer part of a protoplanetary disc. 
The initial conditions are taken directly from a global disc simulation, and a realistic 
equation of state is used to follow the clumps as they contract over several orders 
of magnitude in density, approaching the molecular hydrogen dissociation stage. The 
effects of clump rotation, asymmetries, and radiative cooling are studied. Rotation 
provides support against fast collapse, but non-axisymmetric modes develop and effi- 
ciently transport angular momentum outward, forming a circumplanetary disc. This 
transport helps the clump reach the dynamical collapse phase, resulting from molecular 
hydrogen dissociation, on a thousand-year timescale, which is smaller than timescales 
predicted by some previous spherical ID collapse models. Extrapolation to the thresh- 
old of the runaway hydrogen dissociation indicates that the collapse timescales can be 
shorter than inward migration timescales, suggesting that clumps could survive tidal 
disruption and deliver a proto-gas giant to distances of even a few AU from the central 
star. 

Key words: planet formation - extrasolar planets - protoplanetary dies - gravitational 
instability 



1 INTRODUCTION 



The combined efforts of RV and transit (see Wright et al.|[2012 Borucki et al.pbll I, mircrolensing ( Sumi et al.||2011 l, and 
direct imaging surveys ( Lafreniere et al.||2010| have revealed a wide range of planetary systems, with planets on both very 
short and very wide orbits. These increasingly large sample-size distributions of known planets, although still biased, can be 
used to constrain planet formation theories. Unfortunately, testable predictions by planet formation theories are not fully 
developed, and in some cases, there may be multiple ways to produce a particular planet ( Boley et al.|2010 1 or to induce the 
observed disc structure (see Muto et al.|2012 Hashimoto et al.|2011 l. To complicate matters further, regardless of the planet 
formation mechanism, planet-planet and planet-disc interactions can lead to large scale transport of planets throughout their 
natal discs (see Baruteau et al.|[2011 Michael et al.|20lT |. 

There are currently two main theories (for a review see Armitagc 2010 \ for the formation of massive planets: core accretion 
(hereafter CA) (Mizuno|1980||Pollack et al.|1996| and gravitational instability (hereafter GI) ( |Boss|1997| |Mayer et al.|2002 1. 
From these mechanisms, it is also possible to derive alternative scenarios, such as the tidal stripping theory ( BTley*^r^L[2*0*To 



Nayakshi n|2010[ ). In each paradigm, the formation of a gas giant planet proceeds in a very different way. In the first case the 
grains in the circumstellar disc aggregate until they reach a critical core mass at which there is a run away accretion of gas 
from the disc onto the solid core. The critical core mass for runaway gas accretion can vary due to conditions in the discs 
(Mov shovitz et aT7||2010| |Rafikov||2011[ ); however, the large population of Neptune-mass planets on short periods discovered 
by Kepler suggests that the critical value may be > 10M® . Under the disc fragmentation paradigm, the giant planet is the 
result of a contraction of a massive self-bound clump of gas that forms due to fragmentation of Gl-driven spiral arms. 

The idea of forming giant planets by disc instability goes back to at least iKuiperl (119511). Renewed interest in disc instability 



as a planet-forming mechanism is largely due to the advancement of Boss ( 1997 \ and to recent discoveries of extrasolar planets 
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on wide orbits ( |Bonavita et al.|2011 Nielsen et al.pOlT Quanz et al.|2012 |. The Toomre (1964) stability criterion measures 
the susceptibility of the disc to growing perturbations due to self-gravity; whenever 



2c s 
GXt r 



< 1, 



(1) 



axisymmetric ring distortions will grow in an infinitesimally thin disc. Here E is the disc surface density, c s is the sound 
speed, and tp is the epicyclic period, which will be close to the orbital period in nearly keplerian discs. For three-dimensional 
protoplanetary discs with finite width, numerical simulations have shown that spiral arms will grow whenever Q < 1.7 ( |Durisen| 
eTaLl[2007] . The mass scale for self-gravitating distortions (Toomre mass), which is not necessarily the fragmentation mass, 
can be expressed as a sound horizon 



GMt/cI W Cstp, 



(2) 



which at plausible radii gives massive giant-planet to brown-dwarf masses within a square Toomre wavelength. Simulations 
have shown that when clumps do form, initial fragment masses can be a factor of a few or more smaller than the Toomre 



mass (Mayer et al. 2004 Boley et al. 2010 Rogers & Wadslcy 2012). After spiral structure develops in a disc due to self- 



gravity, spiral arms could regulate the instability and reach a marginally unstable state, with mass redistribution and energy 
dissipation balancing the effects of disc cooling. In the local limit, Gammie (20011 showed that for a cooling time t coo \, 



will lead to an effective |Shakura fc Sunyaev| ( |1973[ ) a viscosity given by 

2 t P 



GIs 



97iT(r - 1) t cool 



(3) 



where T is the (2D) adiabatic index. If GIs cannot stabilize the disc, then spiral arms can fragment into bound clumps. The 
cooling time boundary for fragmentation, which sets a maximum a in the local limit ( Rice et al.|2005 1, is still being debated 
in the literature ( Meru fc BatepOll Lodato & Clarke 2011, Paardckoopcr 2012), but appears to be within a factor of a few 

^cool • 

Now consider the conditions in a disc that are favorable to disc instability. For T = 100 K and tp ~ 12 yr, E > 3600 g 
cm -2 for Q < 1.7. For comparison, a minimum-mass solar nebula has a surface density E ~ 150 gem -2 at 5.2 AU, although 
we also remind the reader that that this value is a minimum and that the assumptions used to derive this mass are suspect 
(for an example, see Desch 2007). For clump formation, higher surface densities than those considered above may be needed, 
as fragmentation seems to require an azimuthally averaged Q < 1.4 (Mayer et al. 20041. The exact threshold will depend 



on the equation of state. Because a low value of Q is needed for fragmentation, it is of interest to know where low values 



will most likely persist. In a Keplerian disc, Q oc r' 



4/2-p-l.S 



for disc radius r. We have assumed that the temperature and 



surface density profiles can be described by power laws with indices q and p, respectively. For a flared disc with q ~ —0.5, 
the surface density must drop more steeply than p ~ —1.75 to prevent the outermost regions of the disc from being the most 
susceptible to GIs. Furthermore, as the surface density drops and the disc becomes colder, the opacity of the disc will also 
drop and cooling can become more efficient, as long the gas has enough emissivity to radiate. Because the local dynamical 
times are long at large radii, cooling times can be much shorter than the local dynamical time, and disc fragmentation, not 
just instability, becomes likely (e.g., |Rafikov|2009| |Boiey|2009| ). 

When clumps do form, they are many orders of magnitude less dense than a giant planet, with radii of a few AU, and 
have significant angular momentum. Although there have been detailed numerical studies of clump evolution, they have been 



limited to spherical quasi-equilibrium models (Helled & Schubert||2009|). Simulations that follow the clump in 3D and within 



the context of a disc are too low resolution to capture the internal dynamics of the clump ( Stamatellos et al.|2007 l. As long 
as clumps remain large and diffuse, it has been suggested that they might lose a significant fraction of their mass and become 
super-Earths as long as they can accrete enough mass in the core ( |Nayakshin||2010| |Boley fc Durisen||2010[ ). However, even 
prior to considering tidal mass loss, the actual mass that goes into forming the planet (or brown dwarf) and the mass that 
contributes to an eventual subdisc (circumplanetary or circum-brown dwarf) is unknown and depends on the efficiency of 
angular momentum transport within the clump. As discussed in more detail below, when the central temperatures of the 
clump reach ~ 2000 K, H2 dissociation will cause rapid collapse of the clump, forming a bound substellar companion (gas 
giants and brown dwarfs). At this point tidal mass loss will be negligible due to a very high central density. 

In the present study we present high resolution 3D simulations of clumps that form through disc fragmentation at 80 AU. 
Their evolution is followed until the contaction approaches hydrogen dissociation and rapid collapse begins. The simulations 
are performed using a new version of Gasoline (Wadsley et al. 20041 that implements an equation of state (herein EOS) 



designed to model the relevant temperature-density regime. In some of the simulations we also include the effect of radiative 
cooling. Section [2] describes our methods. Section [2 . 3| shows simulation results, and discussions and conclusions are presented 
in section [4] A detailed derivation of the EOS is given in the appendix. The purpose of the present work is to follow the 
evolution of clumps in order to build a self-consistent picture of the types of objects that GI forms; this is a rich area of 
research, and this study represents an initial step toward this goal. 



© 0000 RAS, MNRAS 000, 000-000 



Collapse of protoplanetary clumps 3 

Table 1. Simulations parameters for the global simulation presented in Boley et al. 2010 and for the high resolution simulations presented 
in this work. 





Mass [Mj] 


Radius [AU] 


N. particles 


Grav. Softening 


N. neighbours 


Global Simulation 


3.7 


7.0 


1.5 xlO 4 


0.5 AU 


32 


Clump Simulations 


3.7 


7.0 


1.5 xlO 5 


0.045 AU 


32 



2 METHODS 

2.1 EOS and Radiative Cooling 

The internal structure of any single clump will go through a wide range of densities and pressures as the object cools and 
contracts, which can change the relative importance of the translational, rotational, and vibrational modes of molecular 
hydrogen during the clump's evolution. In particular, the onset of H2 dissociation in the clump's core will lead to dynamic 
collapse of the entire structure. For these reasons, we have modified the SPH code Gasoline ( [Wadsley et al.|2004 1 to include 



an EOS for a mixture of atomic and molecular hydrogen in chemical equilibrium (see appendix A). The EOS is taken to be 

= K b Tp(l + a) 

2m P ' V ; 

with P being pressure, kp the Boltzmann constant, T temperature, p density, and a — a(T,p) the dissociation parameter, 
given by the number of protons in the atomic case over the total number of protons, and mp proton mass. In Figure |T]), we 
show the dissociation parameter as a function of density and temperature. Dissociation becomes more likely as the temperature 
increases and as the density decreases. In order to highlight the regime of interest in this context, the evolution of one of our 
simulations is also shown. 

For a clump to evolve from the initial low-density, low-temperature state to the H2 dissociation threshold, it must radiate 
away energy. In the following, we present simulations with and without cooling. The latter is referred to as the adiabatic case, 
and serves as a base for comparisons. We implement cooling in Gasoline for these simulations with a simplified prescription 
that uses local gas conditions only (as done in |Boley et al.|[2~010| ). Let the energy loss per time per volume be given by 

A = (36 7 r) 1 /^(T 4 -T min )-^ T (5) 

where r is the optical depth, s = (m/p) 1,/3 and a is the Stefan-Boltzmann constant. Here, a minimum background temperature 
T mm = 10 K is assumed, which is the reference background radiation field for the ISM. The clump explored here forms in the 
outermost regions of the disc, where the tempertature can become dominated by the radiation background, which is what we 
assume. The optical depth is given by r = pns, where the opacity is approximated as 

(6) 

for T 1 — T if T > T mm and T' = T mm otherwise. While equation (JsJ is only approximate, it allows us to capture the general 
behavior of radiative cooling. Namely, cooling is most efficient at an optical depth r ~ 1. The opacity law used here is also very 
approximate, but is chosen to be monotonically increasing with temperature for two reasons: (1) Only a local optical depth is 
used, where the proper optical depth is integrated along a path. (2) The geometry of the problem places the densest material 
at the highest temperatures. Although sublimation of ices, organics, and dust can all cause the opacity to drop suddenly, the 
integrated optical depth for radiation to leave the clump from any of these regions will likely be large. The simple form for 
the opacity permits fast evaluation of radiative cooling, allowing us to focus on investigating the effects of the EOS on the 
clump's evolution. The effects of radiative transfer in addition to the new EOS will be investigated in future work. 



2.2 Initial conditions 



We present four high-resolution simulations of gaseous clumps that have formed through disc fragmentation. Clump initial 
conditions (ICs) are taken directly from global simulations of a fragmented protoplanetary disc ( |Boley et al.|2010 |, but with 



the resolution increased by a factor of 10 (see below for a description of the procedure). The clump was extracted from the 
global simulation when the fragment's central density was one order of magnitude larger than the background value. Its 
radius was taken to be twice the bound radius. We selected a clump that has a stellar separation of about 80 AU and is not 
distorted due to tidal effects from other close clumps. This choice ensures that the clump has a simple morphology and that 
it is not strongly affected by the presence of other objects or the central potential (see Figure |3|. By using ICs based on the 
results of global simulations, we are confident that the initial conditions of the presented simulations are self-consistent with 
the formation of clumps by disc instability. We refer to the ICs that are produced by this direct extraction method as IC1. 
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As the aim is to follow the collapse of the object through several orders of magnitude in density, the simulations are quite 
computationally expensive. 

We increased the effective resolution in the extracted clump by resampling the "parent" particles (i.e. those taken from 
the global simulation) by ten "child" particles, randomly distributed within a volume defined by the SPH kernel smoothing 
length h determined using 16 neighbors. The child particles inherited equal fractions of the mass from their parents as well 
as their velocity. However, in order to preserve the multiphase structure of the clumps, the hydrodynamic quantities were 
distributed to the child particles using a standard SPH scatter scheme. The temperature has been rescaled to take into account 
the change in the mean moleclar weight fi and in the adiabatic index 7 between the global simulation and the high resolution 
ones used in this study. To simplify the software development, we reused the neighbour-finding methods of the well-known 
group-finding code SKID (Stadcl 2001) and tailored them to our needs. 

The ICs produced by direct extraction can be used to produce a second, quieter set of ICs. This second set (IC2) is 
produced by symmetrizing the mass and the temperature profiles of IC1 and by setting the initial velocity field, which is 
mostly rotational, to zero. This gives us a spherically symmetric IC that can be used to isolate the effects of rotation on 
clump contraction. For each IC, one simulation is run adiabatically (no cooling) and a second is run with the radiative cooling 
described in section 2. The adiabatic case will help to isolate the EOS effects from clump cooling during the early stages of 
the system evolution. 

We check whether the tidal potential of the central star can affect the clump evolution by running an adiabatic simulation 
with IC1 for ~ 10 of its internal dynamical times with the central star explicitely included. We found negligible differences 
between clump evolution with and without the central star. At the time that the clump was extracted from the global 
simulation, the clump is at large stellar separation and well within its Hill sphere. Tidal effects only play a role in the outer, 
low-density layers of the clump, which appear to be unimportant for the evolution of the clump as studied here (see section 
[3] for a detailed analysis). 

The extraction method used to create the initial conditions prevents us from studying the effects of gas accretion from 
the disc onto the clump during its evolution. This is a limitation for these simulations. Gas accretion potentially plays an 
important role in the evolution, although it is not clear in which direction it would lead. As a first approximation, it would 
speed up the collapsing timescale as it increments the clump mass. At the same time, though, this process would influence 
the density and the temperature of the object, leading to a different cooling rate which is not easy to predict due to the 
non-linear nature of cooling. In order to correctly quantify the effects of gas accretion, then, long term simulations on large 
scale are needed, which are beyond the aim of this work. 

Table [1] summarizes the numerical parameters used in this study, including the initial clump mass and radius. Figure [4] 
shows the temperature, density and cumulative angular momentum profiles of the initial clump resampled at higher resolution 
(see next section). The angular momentum barrier Ri, is the radius that the object would have if it were rotationally supported, 
where Rb(< r) = J(< r) 2 /GM(< r) with specific angular momentum inside r J(< r). As the initial radius is 2.5 AU and the 
corresponding angular momentum barrier is 0.17 AU, the clump is only partially rotationally supported. This is confirmed by 
the initial ratio between rotational and gravitational energy, which is 0.18. 

Figure [5] shows the different components of the energy as a function of the radius for the initial condition (kinetic, thermal 
and gravitational). It is evident that the the gravitational component is dominant and that the clump is out of equilibrium. 
The virial equilibrium condition 2(i?i cm + E^ e ) = — -EgraJ^jis not completely fulfilled, as the initial condition is missing 4% 
of the internal energy needed for it to be in virial equilibrium. This out-of-equilibrium condition is due to the low resolution 
of the global disc simulation from which we are extrapolating the clump, as in that condition the clump was not allowed to 
properly collapse because of a gravitational softening of 0.5 AU. This artificial out-of-equilibrium state leads to a fast initial 
collapse phase, which lasts for 4£ c jy n . The clump evolution is then self-consistent after this initial transient. 



2.3 Determining a Resonable Clump Resolution 

The resolution of the clump as taken from the global simulation has been increased by a factor of 10 through particle splitting. 
This value of 10 has been selected based on the results of a resolution study, where IC2 has been evolved adiabatically at 
different resolutions. The initial temperature is kept constant for all particles that are split from their parent. The new particle 
mass is rescaled so that the total mass is constant between different resolution runs and the softening (ft nsm ) is rescaled so 
that ftjjgjjj x (iVpar) 1 '' 3 is constant. Figure [g] shows the evolution of the half-mass radius and of the inner density in different 
runs; there is convergence for both quantities at increasing number of particles. The temperature profiles for three runs at 
different resolutions at different times are shown in Figure[7] In this case, the convergence at high resolution is clear. Moreover, 



1 Note that only the translational component of the thermal energy goes into pressure support. If the total thermal energy is used, then 
the virial condition is dependent on the volume-averaged adiabatic index, such that 2BT,j ne +j c + 3(7 + l)^thermal = — ^grav- 
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in the low-resolution case, the inner part of the clump does not reach as high of temperatures as found in the high-resolution 
runs, as the simulation is not resolved well enough to follow the collapse. 



3 RESULTS 

Using the simulations presented here, we follow the evolution and fate of clumps as they evolve toward second core formation. 
This evolution is strongly correlated to the rotational state of the object, so our analyses first focus on the effects of rotation 
and angular momentum transfer due to the growth of structures in the clump. These processes affect the evolution of the 
clump, which we quantify by studying its half-mass radius, inner density and temperature for the simulated timescale. Lastly, 
we give an estimate of the clump entropy and compare it with the values usually assumed in long-term contraction and cooling 
simulations. 



Radius and angular momentum evolution 

We follow the evolution of the clumps in four different cases: IC1 and IC2 are evolved using the adiabatic and cooling cases. 
The evolution of the half-mass radius (Figure [8]l shows that the initial rotation of the clump (IC2) helps to prevent the 
collapse. This is a trivial result, but note how this effect is quite drastic and fast. In only a few dynamical times the difference 
between the two initial conditions is evident. Naturally, including a cooling term exacerbates the initial collapse. 

The presence of initial asymmetries in the IC1 simulations leads to the development of Fourier modes, which are due to 
low-amplitude spiral structure. They are presumably initially seeded by the self-gravity of the disc and then amplified during 
the collapse phase by the clump dynamics. The stellar potential does not play a role during this phase of the evolution. 

The strengths of the asymmetries are given by the global Fourier amplitudes, 

(a? U/2 
n J poujaujaz 

with 



pcos{mij>)ijjdujdzdcj), and (8) 
b m — I ps\n(m<j>)u)dwdzd4>. (9) 



Integration is extended out to 7 AU. From this analysis it is clear that the clump develops an m-2 mode, which becomes 
weaker with time (see Figures [9) |10| and 111. These structures move the angular momentum outward, as can be seen in 



Figure |12| where, after the fast initial collapse, the cumulative angular momentum profile appears to be flattened in the inner 
part. When the modes start to become unimportant, this angular momentum transport becomes weaker, with the clump 
developing a spherical core (with only 46% of the initial mass collapsed into the core in the cooling case) surrounded by a 
rotating circumplanetary disc. In the adiabatic simulation, although the m-2 mode lasts for a longer time, the effects due to 
the angular momentum transfer are weaker compared to the cooling case. The total angular momentum is conserved in both 
simulations up to 



Figure 13 shows the evolution of the ratio between rotational and gravitational energy T/||W|| for IC1 in the adiabatic 
and cooling cases. In both cases this ratio initially increases, as the collapse is very fast and wins over the angular momentum 
transport. After the first few dynamical times, the collapse becomes slower and the m-2 modes move angular momentum 
outward. Eventually the T/||W|| stops increasing and instead decreases. This phase is faster in the adiabatic case as the 
clump's collapse is less efficient. It is worth pointing out that the initial condition IC1 is not spherical, so the values of T/|| W| 
are not a straighforward indicator of the modes evolution of the clump. Indeed, we detect the presence of several fourier modes 
throughout the collapse, although the ratio hardly reaches 0.274, which is the threashold value for a bar mode in a spheroidal 



( Durisen et al.|1986 l 



While the clump is largely a spheroid, it can also be decomposed into core- like and disc-like structures. For the disc portion, 
the Toomre Q (1964) parameter can be used to explore the susceptibility of the system to the growth of non-axisymmetric 
structure. However, this must be taken with caution, as portions of the disc have significant pressure gradients, making a strict 
application of the Toomre Q difficult. Figure [T4| shows the evolution of Q in the cooling case IC1 outside the half-mass radius 
for five different times. Although early in the clump's evolution the clump's Q value approaches the instability threshold for 



a thick, 3D disc (e.g. Mayer et al. ( 2004 1), the disc evolves toward a stable state. 



At the final stage of the simulation, the protoplanetary disc can be defined as the region between the clump radius and 



2 The conservation of the total angular momentum can not be seen from Figure [l2| as the x-axes are only up to 1.5 AU, so the component 
from the outer part is not shown. 
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Table 2. The last density of the clump as taken directly from the simulations, as well as extrapolations for two evolutionary timescales 
and two thermodynamic values just priori to dissociation and collapse. See the text for further details. 



Quantity IC1: adiabatic IC1: cooling IC2: adiabatic IC2: cooling 

Last density [g/cm 3 ] 4.20 X 1CT 9 1.57 X 1CT 7 4.28 X 1CT 7 4.30 X 10~ 6 

1300 K Timescale [yrs] 2.3 x 10 4 2.2 x 10 3 1.38 x 10 2 

Dissociation Timescale [yrs] 6.3 X 10 3 6.3 X 10 2 

Extr. density [g/cm 3 ] 7.3 X 10~ 6 1.3 X 10~ 4 

Extr. Specific Entropy [k B ] 14.8 15.2 



one third of the Hill radius, which is the maximal disc extent found in previous works (Ayliffe & Bate 2009; Martin & Lubow 



2011 1. Although these previous works address clumps formed via CA, with gas accreting onto solid cores, the disc which forms 



in the simulation herein presented has very similar properties. Figure [15] shows the disc's morphology; the disc gas mass is 
46% of the total mass, giving the clump and disc having comprable mass. Moreover, the clump velocity is sub keplerian in 
the core-like structure due to pressure support, while it approaches keplerian values in the disc-like structure. The ratio H/r 
where H is the disc hight is comparable to the values found in previous works. 



Density and temperature evolution 

To understand if a clump is going to survive inside a disc, it is essential to study the evolution of the clump's inner temperature 
and density, where we use inner to refer to the values that are averaged over the particles inside one gravitational softening 
from the center of the clump. This average is used because the gravitational softening is larger than the particles' smoothing 
length. This also implies that the presented simulations underestimate the contraction, regardless of their high-resolution. 
The evolution of inner density and temperature are shown in Figures [l6| |17| [l8l and |19| The density profile evolution in both 
the adiabatic and cooling cases shows that a fraction of the mass lies outside the Hill radius (13.1 AU in the initial condition), 
which will be stripped away due to the interaction with the host star, reducing the mass of the clump. This process will 
happen on a timescale comparable to the rotational time of the material at the Hill radius, which is (for our IC1) about 760 
yr. This timescale is longer than the duration of the simulations, so neglecting the star is not expected to alter the evolution 
of the high-resolution simulations. However, the timescale is shorter than the contraction timescale, which means tidal effects 
could still play some role in the evolution of the clump over the contraction timescale. This effect is expected to be small for 
the conditions considered here, as the clump becomes highly concentrated over the duration of the simulation. 

The simulations with cooling show that the clump will eventually reach the second core collapse. It is not possible to state 
with certainty the same for the adiabatic simulations, as it is not clear which physical phenomena will lead to the collapse 
after the removal of the m-2 mode. Without cooling, the clump should eventually reach an equilibrium. In the cooling case, it 
is possible to estimate the timescale for reaching the dissociation of molecular hydrogerj^] by extrapolating temperature and 
density (the first quantity is extrapolated linearly, while the latter using a parabolic function). For IC1 the extrapolation is 
done using the last 15 t^y n , while for IC2 it is done using the last 1.5 i^yn' Although it is initially unclear whether such 
extrapolation is reasonable, as after the 3D simulations end, the clumps must still go through a wide range of temperatures 
and densities before collapsing to the second core, the simple extrapolation does place the inner values on a sensible trajectory 
(see Fig. 1). The speed of collapse will be dependent, in part, on the adopted cooling approximation for radiative cooling, but 
also on the opacity and metallicity of the clump, where for a given metallicity the opacity is only known to factors of a few. 
These extrapolations are crude, but we expect the resulting timescales to be within the range of plausible evolution scenarios 
for realistic clumps. See table [2] for the estimated timescales and physical quantities. 



Entropy evolution 

It is usually assumed that CA and GI lead to the formation of protoplanets with different properties. In the first case, called 
the cold model (s/fcg < 10), these objects are supposed to have a lower specific entropy than in the second case, the so called 
hot model (s/fcs > 10) (Spiegel & Burrows 2012). With our simulations we are able to give an estimate of the actual value 
of the specific entropy when the clump reaches R = lORj by extrapolating the entropy evolution (see table |2J. The specific 
entropy of the clump changes very slowly with time, and the extrapolated value agrees with the hot model range found in 
previous works. 

3 In this work we consider the clump to reach the second core collapse when the dissociation parameter is 10% in the inner part of the 
clump. This is a safe assumption as the run away process due to dissociation of molecular hydrogen actually starts when a is only a few 
%. 
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The evolution of specific entropy changes depending whether the inner or the outer part of the clump is considered. It is 
interesting to notice that, in all the performed simulations, the specific entropy of the inner part of the clump decreases with 
time, while it increases for the outer part. This points towards the existence of a redistribution mechanism of the entropy 
between the clump and its envelope. Figure [20] shows the evolution of the specific entropy profile for IC1 in both the adiabatic 
and cooling case. Overall redistibution is evident, although the simulation with cooling develops a more complex specific 
entropy profile with a peak and a trough at r p = 0.48 AU and r± = 2.38 AU, respectively. This effect is due to the more 
effective collapse seen in the simulations with cooling, so that at r = r p the density has significantly decreased compared with 
the adiabatic case (factor of 5), which leads to a local entropy increase. 



4 DISCUSSION AND CONCLUSIONS 

Gravitational instabilities in the outer regions of circumstellar discs can lead to disc fragmentation and the formation of 
clumps, but at the moment it is unclear whether these clumps can survive to become gas giant planets or other bound 
objects. In this work, we have analyzed the initial collapse of realistic clumps using a set of high-resolution, 3D simulations 
that include an EOS that is appropriate for clump densities and temperatures. Our results allow us to estimate the long-term 
evolution of these clumps, which can be used to address their survival and whether they can continue to contract to form a 
companion (in this case, a gas giant planet). 

Our timescale estimates for the collapse are usually smaller than the values found in previous work. In Helled & Schubert 
(2008) the timescale estimated to reach an inner temperature of 1300 K (value at which the dust component evaporates) 
is 1.6 x 10 4 years for a clump with 3Mj. It is possible to make a comparison between this case and our IC1 simulation in 
the cooling case as it implements a similar physics and as the initial conditions have comparable values for both mass and 
luminosity (see Helled k, Bodenheimer|2011 1. Our result for estimated timescale to reach the same inner temperature shows 



that a 3D treatment of the clump evolution describes a faster colla pse (see Table [2^. Th is difference is due to a combination 
of factors, above all the different treatment for contraction. Indeed [Helled fc Schubert] ( |2008| implements a ID quasi-static 
model, which means that the collapse goes through a series of equilibrium states. In this way if reactions take place on a 
dynamical timescale that is smaller than the sound crossing time, the shells are not able to communicate and each of them 
has to wait for the evolution of the others to react to it, while non-axisymmetric and dynamical instabilities can occur in our 
3D simulations. 

The estimated timescales presented in this work have important implications for the evolution and survival of clumps 
formed by disc instability. The clump contraction timescales and central mass concentrations determine whether fragments 
can survive tidal stripping forces as they move throughout the disc. It is possible to estimate the minimum distance from the 
central star that the clump can reach before tidal forces will overwhelm the clump's stability and prevent collapse to a second 
core. To calculate this, use the Hill radius definition: 

«min = ^f' M ^V (10) 



v m clump / 

with Af g ^ ar = O.SMq. Using the results from the more realistic simulation (IC1 with cooling), i.e., the values from the last 
output, we find a mm = 5.97 AU. It is also possible to extrapolate the half-mass radius of the clump to its value as the clump 
reaches collapse to a second core by R = R\ as i (Plast/^extr)*'' 3 — 0-32 AU, so that a mm = 2.52 AU. This means that if the 
clump collapses before getting closer to the star than a mnl , it can become bound enough to survive further tidal effects. This 
is expected to be the case, as the migration timescale for such an object has been found to be of the order of 10 4 yr (compare 



Baruteau et al. 2011 Michael et al. 20111, which is longer than the extrapolated time to molecular hydrogen dissociation 
in the clump. There is the possibility that the clump gets disrupted at pericenter passages in the early stage of migration 
through clump-clump scattering, clump-spiral arm excitation, or birth on an eccentric orbit. If migration remains smooth, 
then the clump can remain well inside its Hill sphere for plausible migration rates. 

Although rapid collapse could be a common evolutionary scenario for clumps, there is still the potential for rich dynamics 
to occur prior to molecular hydrogen dissociation. This is best illustrated by noting the total angular momentum of the 
IC1 clump with cooling at the end of the simulation, which is L — 5.4 x 10 47 g cm 2 s _1 . This result is in agreement with 



the calculations shown in Machida et al. (2008) which give a slightly larger value for L (less than factor of 2) for a clump 



with the same mass. In both cases, the total L is two orders of magnitude larger than the angular momentum estimated 
for Jupiter (Lj — 4.14 x 10 45 g cm 2 s _1 ). Moreover, as the clump evolves, the amplitude of the fourier modes get weaker in 
these simulations, leading to a less efficient removal of the angular momentum from the clump. This implies that there has 
to be a second mechanism, later in the evolution of the clump, in order to match the simulation results for L with Jupiter. 



4 The luminosities are evaluated at half-mass radius in order to exclude the circumplanetary envelope. The initial value is 4.5 X 10 29 erg/s, 
while the end value for the IC1 simulation with cooling is 4.0 X 10 28 erg/s. 
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This mismatch could, for example, be due to the neglect of magnetic effects. The temperature and density regimes that the 
clump will experience do allow for possible development of some magnetic drag, due to thermal ionization, that would lead 
to angular momentum loss (compare with Perna et al.|2010 l. The presence of magnetic fields in the protoplanet could also 
lead to a transfer of its angular momentum to a circumplanetary disc due to coupling of the planetary dipole field lines to 
the disc fluid, as discussed in |Takata fc Stevenson ( 19961. Although the Takata and Stevenson derivation applies only for the 
very last stage of the planet formation, when its radius is only a factor of ten larger than its final value, the mechanism they 
describe can lead to the decreasing of L by a factor of 3 — 4 in this very last phase. Self-gravitating instabilities could also 
become important again as the clump contracts and the core spins up. If T/||W|| once again reaches ~ 0.27, then dynamic 
bar instabilities can be rejuvenated. 

The simulations that have been performed allow us to confirm that the clumps formed via GI are described by the so 
called hot state model. As stressed in |Spiegel fc Burrows| ( |2012| , this is particularly interesting as it can allow the next 
generation of observational surveys to discriminate between the GI and the core accretion model. 

The present work points toward a more central role for the GI theory in the direct formation of gas giant planets, as 
the contraction of disc instability clumps may be very rapid, at least for some if not many conditions. The possibility that 
protoplanetary clumps retain most of their mass even when they reach distances close to the star as a result of inward 
migration opens new scenarios for the origin of close-in extrasolar plantes. Indeed, while previous work has suggested that 



partial stripping and concurrent core formation could turn clumps into super-Earths and/or Neptunes ( |Boley et al. 2010 



Nayakshin||2010| |, our findings suggest that in principle a fraction of clumps formed at large radii by GI can become giant 



planets on close-in orbits. Nevertheless, it is worth pointing out that this result comes from extrapolated values, and not from 
a self-consistent study of a clump's contraction over the full range of its first core phase. This work is a first step towards a 
full description of clump contraction, which appears to be a very complex and computational demanding problem. Note also 
that the EOS herein implemented does not really affect the clump evolution in the presented simulations, but it is expected to 
play a major role in the later evolution. A more realistic cooling prescription using radiative transfer, which is only mimicked 
here, as well as simulations taking into account gas accretion from the disc onto the clump, will need to be included in future 
models. These points will be addressed in future work. 
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Figure 1. Dissociation parameter a as a function of temperature and density. The x and y axes are in log scale (cgs units). The white 
solid curve represents the evolution of the clump from the IC1 simulation in the cooling case, the white dotted curve is the extrapolated 
evolution. 
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Figure 4. Temperature (on the left), density (on the middle) and cumulative angular momentum (on the right) profile of the initial 
clump. The gray dots represents shells containing 10% of the mass. 
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Figure 5. Cumulative energies (gravitational, kinetic and thermal) in log scale (cgs units) as a function of the radius (in AU) for the 
initial condition. 
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Figure 7. Temperature profiles for runs at different resoltions and at different times: from left to right t = lT c i yn , t = 1.5T<j, 
t = 3.0Td yn . As in the previous figure, N is the number of particles in the high-resolution case. 




Figure 8. Evolution (in dynamical times) of the half-mass radius (in unit of its initial value) for IC1 (left) and IC2 (right) in adiabatic 
(solid line) and cooling (dashed line) cases. 



© 0000 RAS, MNRAS 000, 000-000 



Collapse of protoplanetary clumps 13 




Figure 9. Fourier analysis of the evolution of IC1 at different times: adiabatic case. The maps show the difference between the density 
and the mean density calculated from the surface density profile of the clump. Axes are in AU, sec the bottom left corner. On the bottom 
right, corresponding amplitude vs Fourier modes 




Figure 10. Fourier analysis of the evolution of IC1 at different times: cooling case. The maps show the difference between the density 
and the mean density calculated from the surface density profile of the clump. Axes are in AU, see the bottom left corner. On the bottom 
right, corresponding amplitude vs Fourier modes 
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Figure 11. Evolution (in dynamical times) of the amplitude of the first four fourier modes in the adiabatic (left) and cooling (right) 
cases using IC1. The solid line corresponds to m=l, the dottcd-dashed line to m=2, the dashed line to m=3 and the dotted line to m=4. 
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Figure 12. Evolution of the cumulative angular momentum profile in the adiabatic (left) and cooling (right) case for IC1. The transient 
phase (earlier than 4t c jy n ) is not shown. 
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Figure 13. Evolution of the ratio between rotational and gravitational energy T/||W|| for IC1 in the adiabatic (on the left) and cooling 
(on the right) cases. 



o . 


t = 0td 

t = 8.9td 






51 11 








t = 17.8 td 






1 1 1 1 

1; 

'1 1 1! 
ii 1 






CO 


t = 53.5td 

t = 89.1 td 










CO 

o 

"* 


A 
/\ 

/ \ 

/ \ 
1 \ 

A * 




1 

II 
it 

1 

' i 

I 

1 !•'! 


li 

!! f 

1 

s!l 1 






CM 


Vn < 


1 Ai 

tr 








o 






-0.5 0.0 


0!5 

Log 10 (R A u) 


1.0 


1:5 



Figure 14. Toomre parameter Q as a function of radius (in log scale) at five different times for the cooling evolution of IC1. The 
horizontal line is the thrcashold 0.7 for a thick, 3D disc (see, e.g., Mayer et a. 2004). 
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Figure 15. Circumplanctary disc properties at the last stage of IC1 simulation in cooling case. The left panel shows the ratio between 
the disc hight and the disc radius, the central panel the gas circular speed to the Keplerian speed, and the right panel the mass enclosed 
within a given radius. All profiles are shown versus clump radius. In each panel, the vertical lines represent the clump core and the disc 
boundaries. 




Figure 16. Evolution (in dynamical times) of the inner density (log scale of cgs units) for IC1 (left) and IC2 (right) in adiabatic (solid 
line) and cooling (dashed line) cases. 
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Figure 17. Evolution (in dynamical times) of the inner temperature for IC1 (left) and IC2 (right) in adiabatic (solid line) and cooling 
(dashed line) cases. 
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Figure 18. Density profile evolution for IC1 in the adiabatic (on the left) and cooling (on the right) cases. 
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Figure 19. Temperature profile evolution for IC1 in the adiabatic (left) and cooling (right) cases. 
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APPENDIX A: EQUATION OF STATE 

In this Appendix we derive the equation of state for a mixture of atomic and diatomic hydrogen. We define the dissociation 
fraction a as 

Nvt 

a=^ (Al) 

with TVh number of atoms and N p total number of protons. 

The atomic and molecular parts are assumed to be in chemical equilibrium. Within the molecular part, however, the two 
nuclear-spin configurations (ortho and para hydrogen) are not in equilibrium, but taken to be frozen at ortho : para = 3:1, 
as there is no efficient mechanism to convert them into each other in dusty clumps (Boley et al.|2007). 



Atomic component 

We now derive the thermodynamic functions for the atomic component. The partition function will be made up of only the 
translational component q t per particle 



where V is the volume and A a t = h/y/ (2TvmHk B T), with h being Planck's constant, tuh the atomic mass of hydrogen, k B 
Boltzmann's constant and T the temperature. From that we can calculate logQii: 

log Q H = log(gf H ) ~ TVh log = Afa log — 13- , (A3) 

where uh is the number density. Now we can apply the definitions for the Helmholtz energy Ah , internal energy Eh , pressure 
Ph and entropy Sh to derive 

A H = -k B T log Qh = -k B TNn log (A4) 
Eh = k B T 2 ^°^ =%NhT (A5) 

° N„,V 2 



p u rr dlog Qh 



_ NHTk B _ k B pT 

— 77 — (, A "J 

N H ,T V m H 



Sh = ^5 = ^ fl (| + Iog « ,. 



Diatomic component 

In the diatomic case, the partition function has a rotational component q r and a vibrational component q v . 

It = Z- (A8) 

'Vol 

qr = gy 4 xg 3/4 =4 /4 4 /4 (e XP (2e r /T)) 3 / 4 (A9) 
q v = l/(l-exp(-0„/T)), (A10) 

where A mo i = h/ y/ {4,-KiriHk B T), B r = 85 K is the critical rotational temperature and Q v — 5987K critical vibrational 
temperature, and 

S E = Y, & + X ) eX P ("JO' + l)©r/T) (All) 
j even 

So = Y, 3(2j + 1) exp + l)@ r /T) . (A12) 

j odd 

Prom this we can derive log Qu 2 '• 

~ 1 log (S E S ) + log « log (1 - exp (-0„/T)) + | % (A13) 

Nh 2 4 A^ lol nH 2 2 T 
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Ah 2 
Eh 2 
Ph 2 

where 
Se 



We can now derive the thermodynamical functions for diatomic hydrogen: 



= -k B TNu 2 
= k B Nu 2 T [>i 



- log (S E (So) 3 ) + § ^ + log - log (1 - cxp (-0„/T)) 



e, 



cMogQH^ jv H2 rfc B 

KflJ ^77 |iV H o. T = 77 = "S-t ™H 2 



©i, 1 O r Se 30 r 5o 

~T expQ v /T- 1 + 4T5^ + 



N H2 k B 



dV 
3 
2 



<,l/4 <,3/4 



+ 10g (1 - e-<-W^)n H2 A^ ol + AT \S E So) T e«W T - 1 



= E 3{j + l)(2j + 1) cxp + l)6 r /T) 

j even 



So = 3 3(J + l)(2j + 1) cxp + \)® r /T). 

j odd 



Dissociation fraction from chemical equilibrium 

We consider our system in chemical equilibrium, so that 

2H H 2 . 

The value of the dissociation parameter a is determined by the equilibrium condition 

[H] 2 _ Q|(int) AL, 



(A14) 
(A15) 



[H2] A« t Qn 2 (int) 

The left side of the equation corresponds to the ratio between the concentration of atomic / diatomic hydrogen: 

[H] 2 _ A>1 V N 2 a 2 IV 2a 2 N p 2a 2 n 

[H2] ~ V 2 ]Vh7 ~~ V 2 (l-a)N p ~ 1/(1 -a) ~~ (1 - a) ' 
where n is the total number density. For the right side of the equation, we recall the functions used above, taking into account 
that <2i(int) is the total internal partition function, which does not consider the translational term and instead has an extra 
term 



(A16) 



e X p(-E g /{k B Tj) 

with Eg ground state energy. Putting everything together, we find that the equilibrium condition is 



2a 2 



F(T) 



with 



F(T) = 



(1 — a) n 
V(^nnk B T) \ 3 W(T)) 1 - e ( - eWT) 



h J 4 /4 4 /4 exp(3e r /(2T))' 
where _Edi s is the dissociation energy. The dissociation fraction as a function of temperature and number density is 

o(T , n) l(_^ + « + 8 ^ 
4 \ n V n 2 n 



(A17) 
(A18) 

(A19) 
(A20) 



The equation of state 

As we know that the pressure is additive, we can find the total pressure P as 
/' /V: • /'h, ----- + 



N H Tk B N H2 Tk B _ k B T (l-a)N p _ k B Tn{l + a) 



V 



V 



~y-aN p - 



(A21) 
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